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Abstract:  This  report  documents  a  Fortran  version  of  a  Pitzer  equations  for  activity  coefficients  and  water  activ- 

chemical  thermodynamic  model  for  aqueous  electro-  ity  calculation.  The  model  includes  both  the  freezing 

lyte  solutions  at  subzero  temperatures,  FREZCHEM2,  (melting)  reaction  pathway  at  fixed  water  amount  and 

which  is  a  further  development  of  the  FREZCHEM  the  evaporation  (dilution)  pathway  at  fixed  temperature, 

model.  The  model  uses  thermodynamic  data  of  Spen-  The  FREZCHEM2  model  can  be  extended  with  respect  to 

cer-Mpller-Weare  that  permit  the  calculation  of  chem-  independent  components,  electrolyte  species,  and  sol- 

ical  equilibria  in  the  Na-K-Ca-Mg-Cl-S04-H20  sys-  ids,  and  if  corresponding  thermodynamic  data  are 

tern  between  -60  and  25°C  at  atmospheric  pressure.  It  available,  the  model  may  be  used  to  compute  chemical 

applies  the  Gibbs  energy  minimization  method  for  equilibria  in  any  systems  that  include  aqueous-solution 

chemical  equilibrium  computation  combined  with  and/or  one-component  solid  phases. 
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FREZCHEM2 

A  Chemical  Thermodynamic  Model  for 
Electrolyte  Solutions  at  Subzero  Temperatures 

MIKHAIL  V.  MIRONENKO,  STEVEN  A.  GRANT,  GILES  M.  MARION,  AND  RONALD  E.  FARREN 


INTRODUCTION 

The  FREZCHEM  model  was  developed  by  Mar¬ 
ion  and  Grant  (1994)  to  calculate  chemical  equilib¬ 
ria  among  aqueous  electrolyte  solutions,  ice,  and 
salts.  The  model  applies  the  Pitzer  equations  for 
calculation  of  aqueous  species  and  water  activi¬ 
ties.  To  find  chemical  equilibrium,  this  program 
solves  sequentially  a  set  of  nonlinear  equations 
that  includes  both  solid-phase  deposition  and  ion- 
pair  formation  using  an  individual  subroutine  for 
every  reaction.  FREZCHEM  uses  data  on  constants 
of  chemical  reactions  and  Pitzer  equation  parame¬ 
ters  published  by  Spencer  et  al.  (1990).  The  results 
of  modeling  show  good  agreement  both  with 
experimental  data  and  with  the  results  of  the 
Spencer-Moller-Weare  model.  However,  the 
FREZCHEM  model  has  some  limitations.  One  is 
convergence  problems  at  high  ionic  strengths  (>15 
molal)  and  at  junctions,  where  new  phases  begin 
to  precipitate.  Another  is  that  addition  of  any  new 
substance  into  this  model  requires  changes  not 
only  in  data  but  also  in  the  program  code. 

The  objective  of  this  report  is  the  further  devel¬ 
opment  of  the  chemical  thermodynamic  model 
FREZCHEM  to  make  it  more  reliable,  universal, 
and  flexible.  The  point  calculation  reliability  was 
improved  by  applying  the  Gibbs  energy  minimi¬ 
zation  approach  to  computing  equilibrium.  The 
thermodynamic  information  needed  for  com¬ 
putations  is  separated  from  the  calculating  rou¬ 
tines.  That  allows  components  to  be  added  to  the 
system  without  code  changes  and  the  program 
code  to  be  applied  for  other  chemical  systems.  It 
should  be  noted  that  the  Pitzer  approach  describes 
most  interactions  in  aqueous  solution  as  electro¬ 


static  and  only  explicitly  recognizes  a  few  chemi¬ 
cal  interactions,  such  as  ion-pair  formation.  This 
is  why  the  system  under  consideration  is  very 
simple  from  the  viewpoint  of  chemical  interac¬ 
tions,  but  it  is  very  complex  from  the  viewpoint  of 
the  influence  of  activity  coefficients  on  the  behav¬ 
ior  of  the  Gibbs  energy  function. 

MATHEMATICAL  ALGORITHM 

The  system  under  investigation  consists  of  the 
following  components:  1)  solid  salts  of  fixed 
chemical  composition  and  pure  ice  (so-called 
one-component  phases),  and  2)  aqueous  solu¬ 
tions  consisting  of  water  and  dissolved  electro¬ 
lytes.  The  applied  algorithms  of  chemical  equilib¬ 
ria  computation  will  be  described  in  terms  of 
these  components. 

The  equilibrium  composition  of  the  system  at 
constant  T,  P,  and  specified  bulk  composition 
may  be  found  by  minimizing  the  Gibbs  energy 
function  of  the  system  under  balance  restrictions. 

Local  minimum  computation 

Local  minimum  is  considered  as  an  equilibri¬ 
um  composition  of  the  system,  in  which  all  exist¬ 
ing  phases  are  specified  before  computation.  The 
Gibbs  energy  function  of  the  system  that  contains 
M  solids  and  aqueous  solution  (water  and  /  spe¬ 
cies)  is  as  follows: 

G  M  o  3 

S  =  -^=  +  Y.V-jnj 

K1  k= 1  7=1 

where  G  =  free  energy  of  the  system, 

n  =  the  molal  quantity  of  components. 


(ij  =  the  chemical  potential  of  species  j, 
jLtvv  =  the  chemical  potential  of  water, 

=  the  standard  chemical  potential  of  a 
one-component  solid-phase  k,  and 
R  =  the  universal  gas  constant. 

The  chemical  potential  of  aqueous  solution  spe¬ 
cies  /  in  terms  of  molality  is  defined  by 

(ij  =  (ij3  +  lnflj  =tf  +  ln(mj  Yj) 

where  p°  is  the  standard  chemical  potential, 

n\ 

m\  =  — •  55.5 1 
"w 

is  the  moles  of  the  ;th  species  per  1  kg  of  water 
(molality),  and  y  is  an  activity  coefficient. 

The  chemical  potential  of  water  may  be  writ¬ 
ten  as 


M-w  =iC  +  ln  av 


where  water  activity  tfw,  according  to  Pitzer 
(1987),  is  defined  through  the  osmotic  coefficient 
of  the  solution  (j)  and  molalities  of  species  by 
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where  W  is  the  molecular  weight  of  water 
(18.0153). 

Accordingly,  the  free  energy  function  of  the 
system  is  as  follows: 

M  n  n  Zn; 
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k-\ 


J 

+  i>j 

;'= i 


Ms  55.51  "l 

p?  +  In 

— - Yj 

J 

l  "w  Vj 

(1) 


ponent  i  in  the  system.  For  the  electroneutrality 
equation  b{  =  0  and  Vy  =  Zj ,  where  Zj  is  the  charge 
of  the/th  component.  In  matrix  notation,  eq  2  may 
be  written  as 

N  nT  =  b 

where  N  is  the  stoichiometric  matrix,  n  is  the  vec¬ 
tor  of  numbers  of  moles  of  species,  and  b  is  the 
vector  of  bulk  chemical  composition  of  the  sys¬ 
tem. 

It  is  convenient  to  solve  the  system  of  linear 
equations  (eq  2)  with  respect  to  P  components, 
including  M(M<P)  solids. 
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(2b) 

and  in  this  way  to  switch  to  new  independent 
components.  In  such  a  manner  the  stoichiometry 
of  other  /  +  1  -  P  components  (vectors  <3j )  and  the 
matter  balance  B  are  now  defined  through  these 
independent  components.  This  operation  allows 
the  number  of  active  constraints  to  be  reduced  up 
to  P-M.  The  thermodynamic  meaning  of  this  lies 
in  the  fact  that  the  chemical  potential  of  a  one- 
component  solid  phase  is  equal  to  the  standard 
Gibbs  energy  of  formation  and  does  not  depend 
on  its  amount,  until  this  phase  is  present.  This  is 
why  the  system  can  be  considered  to  be  open 
with  respect  to  this  component. 

It  is  obvious  that 

nw  >  0  and  nj  >  0.  (3) 

Minimization  of  the  function  in  eq  1  under  the 
constraints  of  eq  2  and  eq  3  can  be  replaced  by  a 
search  of  the  extremum  of  the  Lagrangian  func¬ 
tion,  which  may  be  written  as 


Mass  balance  constraints,  including  the  electro¬ 
neutrality  equation  if  necessary,  may  be  written 
as  a  system  of  linear  equations: 


M 


7+1 


®(n,  A,)  =  g(n)  +  (Bk  -  X  «kj«j ) 


k= 1 


M+l+J 
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M 

where  P  is  the  number  of  independent  chemical 
components  in  the  system,  and  Vy  is  the  number 
of  moles  (stoichiometric  units)  of  independent 
component  i  in  one  mole  of  component/.  b\  repre¬ 
sents  the  number  of  moles  of  independent  com¬ 
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where  X  is  a  Lagrangian  multiplier.  It  can  be 
shown  (Karpov  et  al.  1976)  that  X  is  the  chemical 
potential  of  the  corresponding  independent  com¬ 
ponent  of  the  system.  In  particular  for  solids, 

=Fk* 
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The  conditions  of  extremum  of  the  Lagrangian 
function  are  found  where  all  first  partial  deriv¬ 
atives  with  respect  to  components  and  to 
Lagrangian  multipliers  are  equal  to  zero.  This 
gives 

30  n  n\ 

=  (nf  +  In  55.51) +ln(-LYj) 

OWj  J  ttw  ; 
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/=l 

There  are  different  approaches  to  searching  the 
Lagrangian  function  extremum.  One  of  them  is 
the  algorithm  developed  by  White  (1958,  1967) 
for  homogeneous  gas  systems,  which  has  been 
further  developed  for  heterogeneous  multiphase 
systems  by  Karpov  (1976).  This  algorithm  was 
applied  by  the  senior  author  of  this  report  for 
computation  of  a  wide  range  of  chemical  equilib¬ 
ria  (Mironenko  1991,  1992)  and  is  build  into  the 
DiaNIK  system  (Khodakovsky  1992).  The  idea  of 
the  method,  as  applied  to  the  system  under  con¬ 
sideration,  is  as  follows.  Equation  4a  may  be 
solved  with  respect  to  ny 


by  Pitzer's  routine  at  every  iteration.  Because  of 
this  it  was  very  difficult  to  reach  the  required  pre¬ 
cision  of  solution  (0.1%),  even  when  special  steps 
were  undertaken. 

Another  approach  is  to  solve  the  whole  system 
ofP-M  +  J  +  1  equations  (eq  4  a,b,c)  iteratively  by 
Newton's  method  for  ny  nW/  and  X^.  This  algo¬ 
rithm  has  been  described  in  detail  by  Harvie  et  al. 
(1987).  It  has  been  successfully  applied  by  Spen¬ 
cer  et  al.  (1990)  for  strong  electrolyte  solution 
modeling,  but  a  working  version  of  the  program 
has  not  been  published.  This  algorithm  also  was 
applied  by  Mironenko  (1983)  for  modeling  fluid- 
rock  interactions  during  hydrothermal  uranium 
ore  formation. 

The  second  partial  derivatives  of  the  Lagran¬ 
gian  function  are  equal  to 

92Q  _—,ifj  =  jl 
dnjdnfl  [{ifjtjl 
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Substitution  of  these  terms  into  eq  4b  and  c  gives 
the  system  of  P-M  +1  equations^ which  may  be 
solved  by  Newton's  method  for  X  and  nw  .  The 
advantages  of  this  approach  are  a  fast  rate  of 
computation,  due  to  a  small  amount  of  variables 
(their  amount  does  not  depend  on  the  number  of 
species),  and  the  lack  of  necessity  to  undertake 
special  steps  to  calculate  species  at  very  low  con¬ 
centrations  or  to  correct  negative  values  of  mass 
for  species  during  iterations.  Unfortunately, 
attempts  to  apply  this  approach  to  brine  systems 
in  combination  with  Pitzer's  routine  have  dem¬ 
onstrated  that  the  algorithm  is  not  tolerant  of  os¬ 
cillations  of  activity  coefficient  values,  provided 


3A,j3A,ji 

This  matrix  of  second  partial  derivatives  is 
known  as  the  Hessian  matrix. 

Activity  coefficients  of  species  and  the  osmotic 
coefficient  are  calculated  at  every  iteration  using 
Pitzer's  model  (Pitzer  1987).  In  FREZCHEM2,  the 
Pitzer  routines  published  by  Marion  and  Grant 
(1994)  were  used,  with  insignificant  changes 
dealing  mainly  with  the  interface  with  data  files. 

The  molal  amounts  of  solids  are  calculated 
after  a  local  equilibrium  has  been  achieved  using 
eq  2a. 
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Search  for  the  equilibrium  phase  assemblage 

Solids 

If  a  molal  quantity  of  solid  calculated  using  eq 
2a  was  negative,  this  solid  phase  was  considered 
as  completely  dissolved,  and  a  new  local  equilib¬ 
rium  without  this  solid  was  computed.  Then  a 
search  of  new  solid  phases  from  the  list  of  possi¬ 
ble  solids  to  be  included  into  the  system  is  under¬ 
taken.  The  criterion  for  the  inclusion  of  phase  K  is 
as  follows: 

M  P 

X  M-k«kK  +  X : Mik  -  <  0.  (6) 

Jfc=l  i=M+ 1 

The  thermodynamic  meaning  of  this  expression 
is  that  the  free  energy  of  chemical  reaction  of  a 
given  solid  substance  formed  from  independent 
components  of  the  system  is  negative,  and  there¬ 
fore  this  solid  is  thermodynamically  stable.  If  the 
condition  (eq  6)  asserts,  this  solid  replaces  one  of 
the  independent  components  of  the  system  with 
which  it  is  linearly  dependent.  Then  the  system  of 
linear  equations  (eq  2a,b)  is  solved  with  respect  to 
this  new  independent  component.  By  this  means 
the  chemical  composition  of  all  components  of 
the  system  will  be  expressed  in  terms  of  this  and 
other  independent  components.  This  procedure 
is  largely  achieved  by  applying  the  Simplex  rou¬ 
tine.  (Simplex  is  a  classic  finite  iteration  method 
of  linear  programming  [Korn  and  Korn  1963].) 
Addition  of  each  solid  phase  reduces  the  number 
of  active  linear  restrictions  by  one.  Calculations 
are  continued  until,  in  the  list  of  possible  phases, 
there  is  no  phase  that  meets  the  condition  in  eq  6. 

Aqueous  solution 

Aqueous  solution  is  considered  absent  in  the 
system  when  the  number  of  active  balance  re¬ 
strictions  (P  -  M )  is  less  than  or  equal  to  one  and 
the  amount  of  water  is  less  than  0.001  moles. 

Special  steps 

Usually,  the  approximate  phase  composition 
of  a  system  may  be  determined  at  the  first  steps  of 
calculation  using  the  Simplex  routine.  Then  the 
exact  equilibrium  composition  may  be  computed 
using  Newton's  method  for  local  equilibria  deter¬ 
mination  and  Simplex  methods  for  addition  or 
substitution  of  solids.  Due  to  the  very  high  non¬ 
ideality  of  brines,  this  technique  collapsed,  and 
some  changes  in  the  logical  pattern  of  calculation 
were  made: 

1.  At  first,  the  system  is  considered  homo¬ 
geneous  (no  solids),  then  the  solid-phase  assem¬ 


blage  is  calculated,  not  simultaneously  during 
one  application  of  the  Simplex  routine,  but  se¬ 
quentially.  Another  phase  is  added  after  local 
equilibrium  with  previously  added  phases  is 
achieved. 

2.  After  the  appearance  of  a  new  solid  phase 
and  before  applying  Newton's  method,  the  cur¬ 
rent  species  concentrations  have  to  be  recalcu¬ 
lated  to  be  in  better  agreement  with  values  of 
independent  component  chemical  potentials.  The 
relation  between  concentration  of  species  and 
values  of  chemical  potentials  of  independent 
components  is  expressed  by  eq  5  and  can  be  also 
treated  in  terms  of  the  free  energy  of  the  chemical 
reaction  of  species  formation  from  independent 
components  of  the  system. 

3.  Because  of  particularities  of  Pitzer's  model, 
to  prevent  wide  fluctuations  during  solution  of 
the  system  of  equations  4a,b,c  by  Newton's  meth¬ 
od,  we  have  to  smooth  changes  of  activity  coeffi¬ 
cient  and  osmotic  coefficient  values,  which  are 
calculated  at  each  iteration  by  Pitzer's  routine. 
We  use  average  values  obtained  at  the  current 
and  previous  iterations. 

4.  At  every  iteration  a  new  approximation  to 
the  solution  is  provided  by  inversion  of  the  Hes¬ 
sian  matrix:  yp+1^  =  yW  +  A*  /  i; .  For  a  homoge¬ 
neous  system  the  value  of  £  is  equal  to  1  and  it 
increases  by  0.5  with  every  new  solid  that  precip¬ 
itates. 

FREZCHEM2  PROGRAM 

A  listing  of  the  FREZCHEM2  Fortran  program 
is  in  Appendix  A.  FREZCHEM2  consists  of  a  main 
program  called  READWRITE  and  seven  subrou¬ 
tines. 

The  READWRITE  program  reads  input  data 
from  the  file  INPUT,  according  to  these  data 
forms  independent  components  of  the  system, 
and  reads  the  temperature  interval  and  tempera¬ 
ture  step  for  freezing,  or  the  water  content  inter¬ 
val  and  water  decrement  at  a  given  temperature 
for  the  evaporation  scenario.  It  calculates  chemi¬ 
cal  potentials  of  the  components  as  functions  of 
temperature,  calls  various  subroutines,  and 
writes  results  of  the  chemical  equilibria  computa¬ 
tion  into  the  file  RESULT. 

Subroutine  CHOICE  is  called  from  the  main 
program  and  chooses  components  that  may  be 
formed  in  the  system  of  given  chemical  composi¬ 
tion  as  well  as  their  stoichiometry.  A  data  file  for 
this  routine  is  the  DATABASE  file. 

Subroutine  SOL  is  called  from  READWRITE 
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and  is  the  main  calculating  routine.  It  computes 
the  equilibrium  composition  of  the  system  at  giv¬ 
en  T  and  specified  mass  balance  by  searching  for 
the  global  extremum  of  the  Lagrangian  function. 
It  forms  the  Hessian  matrix  and  calls  various  sub¬ 
routines.  Results  of  computations  return  to  the 
READWRITE  program  to  be  written. 

Subroutine  SIMPL  is  called  from  the  SOL  sub¬ 
routine  to  enter  new  solid  phases  into  the  system 
and  to  invert  the  stoichiometric  matrix.  This  calls 
the  GG1  routine,  which  is  the  short  version  of  the 
Gauss-Jordan  matrix  method. 

Subroutine  PITZER  is  called  from  the  main 
program  to  choose  data  for  Pitzer  parameters  and 
to  calculate  them  at  various  temperatures  if  the 
freezing  scenario  has  been  chosen.  It  is  also  called 
from  the  SOL  routine  to  calculate  activity  coeffi¬ 
cients  of  species  and  water  activity,  using  the 
Pitzer  model  at  every  iteration  while  minimizing 
the  free  energy  function  by  Newton's  method.  In 
FREZCHEM2,  the  PITZER  subroutine,  as  well  as 
the  INTERACT  subroutine  published  by  Marion 
and  Grant  (1994),  were  used  with  only  insig¬ 
nificant  changes  dealing  mainly  with  the  inter¬ 
face  with  data  files. 

Subroutine  INTERACT  calculates  the  higher- 
order  electrostatic  interactions  for  the  Pitzer 
equations. 


Figure  1.  Flowchart  of  the  FREZCHEM2  model. 


Subroutine  GG  is  called  from  the  SOL  subrou¬ 
tine  and  solves  the  system  of  linear  equations  by 
the  Gauss-Jordan  method. 

The  principal  flowchart  of  FREZCHEM2  is 
shown  in  Figure  1.  The  FREZCHEM2  model  is  a 
universal  model  that  may  be  used  to  compute 
chemical  equilibrium  in  any  system  consisting  of 
one-component  solids  and/or  aqueous  solution. 
For  this  goal,  only  additions  in  the  data  files  are 
needed. 

Data  files 

Files  that  contain  the  information  needed  for 
calculations  and  some  remarks  are  listed  in  Ap¬ 
pendix  B. 

File  DATABASE  contains  a  list  of  independent 
components,  which  can  be  taken  into  account, 
and  lists  of  aqueous  solution  species  (cations,  an¬ 
ions,  neutral  species  including  water)  and  solid 
phases,  consisting  of  given  independent  compo¬ 
nents  as  well  as  their  stoichiometry.  For  the  con¬ 
venience  of  users,  the  same  species  numbers 
(coding)  were  used  as  in  the  FREZCHEM  model. 
If  necessary,  additional  information  may  be  add¬ 
ed  for  independent  components  as  well  as  for 
species  and  solids.  Requested  formats  for  enter¬ 
ing  new  data  could  be  taken  from  the  listing  of 
the  READWRITE  program  (Appendix  A). 

File  TABLE1  represents  Table  1  of  Spencer  et  al. 
(1990),  which  includes  constants  for  the  Debye- 
Hiickel  model  parameter  and  for  the  binary 
interaction  parameters  as  a  function  of  tempera¬ 
ture  (K).  File  TABLE2  represents  Table  2  of  that 
paper  for  mixed-salt  parameters.  File  TABLE3 
contains  coefficients  for  calculation  of  free  ener¬ 
gies  of  chemical  reactions  of  formation  for  solids 
and  ion  pairs  from  aqueous  solution  species  and 
liquid  water  as  a  function  of  temperature,  using 
equations  of  the  form  published  by  Spencer  et  al. 
(1990): 


-AG 

RT 


=  d\  +  afT  +  afF^  +  afT^ 


+  fl3/T  +  a4ln(T)  (7) 

and  represents  a  copy  of  Table  3  from  their  paper. 
In  this  convention,  free  energies  of  cations,  an¬ 
ions,  and  liquid  water  are  taken  to  be  equal  to 
zero  at  any  temperature. 

Program  input  and  output 

Input  to  FREZCHEM2  is  through  the  file  IN¬ 
PUT,  which  contains  the  molal  amounts  of  inde- 
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Table  1.  FREZCHEM2  model  input  for  freez¬ 
ing  seawater  from  0°C  down  to  -40°C  with  a 
2°C  decrement. 


SMW 

seawater 

Title  of  the  task 

0.48695 

Sodium  (mol /kg) 

0.01063 

Potassium  (mol /kg) 

0.00953 

Calcium  (mol/kg) 

0.05516 

Magnesium  (mol /kg) 

0.56818 

Chloride  (mol /kg) 

0.02939 

Sulfate  (mol /kg) 

0.0 

Carbonate  (mol /kg) 

0.0 

Hydrogen  (mol/kg) 

273.15 

initial  temperature 

1 

freezing  (2  for  evaporation) 

233.15 

final  temperature  (final  amount  of  water  for 
evaporation ) 

2.0 

temperature  decrement  (water  decrement  for 
evaporation) 

pendent  components  (presently  Na+,  K+,  Ca+2, 
Mg+2,  Cl-  ,  SOJ2/  and  CO32  are  included)  per  1 
kg  of  water.  Thus,  an  initial  amount  of  water  in  the 
system  is  equal  to  1  kg  or  55.51  moles.  If  the  molal 
quantity  of  an  independent  component  is  equal 
to  zero,  all  substances  in  the  file  DATABASE  that 
contain  it  will  be  ignored.  The  program  calculates 
a  charge  balance  and,  if  it  is  different  from  zero, 
proposes  to  add  some  amount  of  any  cation  or 
anion  depending  on  the  value  and  the  sign  of  the 
charge  balance,  after  which  the  program  starts  its 
work.  FREZCHEM2  is  able  to  calculate  both  1)  a 
cooling/heating  scenario  (in  this  case  one  enters 
initial  and  final  temperatures  and  temperature 
decrement/ increment)  and  2)  an  evaporation/ 
dilution  scenario  at  constant  temperature  (in  this 
case  it  needs  a  temperature,  a  water  amount 
decrement/ increment,  and  a  final  mass  of  water). 
The  model  is  also  able  to  calculate  the  process  of 
ice  evaporation.  An  example  of  input  (file  INPUT) 
is  present  in  Table  1  for  seawater  freezing  from 
0  to  -40°C  for  a  2°C  temperature  decrement.  The 
explanations  are  written  in  italics. 

Output  from  the  program  is  to  the  RESULT  file. 
As  examples  of  output,  results  of  computation  of 
seawater  freezing  at  228.15  K  (-45 °C)  and  at 
218.15  K  (-55°C)  are  given  in  Table  2.  Output  of 
evaporation  of  seawater  at  273.15  K  (CPC)  down 
to  50  g  of  water  is  given  in  Table  3. 

Distribution  of  the  independent  components 
among  solids  and  solution  phases  is  given  in  the 
BALANCE  table  at  the  bottom  of  the  output  table. 
The  last  column  of  this  table  has  been  printed  to 
show  the  equivalence  of  input  and  computed  bal¬ 
ances. 


Verification  of  the  model 

To  verify  the  program,  phase  diagrams  from 
Spencer  et  al.  (1990)  and  point  computations  from 
Marion  and  Grant  (1994)  were  recalculated.  The 
model  reproduces  these  computations  with  good 
accuracy.  Table  4  shows  the  temperatures  at  the 
appearance  of  solids  during  seawater  freezing, 
taken  from  Spencer  et  al.  (1990)  with  an  added 
column  obtained  by  the  FREZCHEM2  model 
using  their  thermodynamic  data. 

It  is  interesting  to  note  that,  according  to  the 
free  energies  of  chemical  reactions  in  the  model,  a 
solid  reaction 

NaCl  •  2H20(cr)  -»  NaCl(cr)  +  2H20(cr  j) 

takes  place  at  temperatures  lower  then  -57.15°C. 
To  verify  this  independently,  the  heat  capacity 
equations  for  these  phases  at  low  temperatures 
are  needed. 
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Table  2.  FREZCHEM2  model  output  for  freezing  seawater  at  -45°C  and 
at  -55°C. 


Seawater  freezing 
temperature  -A5.00°C  (228.15  K) 

SOLID  PHASES 


N 

Phase 

Moles 

-G/RT 

1 

53.52785 

-0.4278 

2 

NaCl*2H20(cr) 

0.42624 

1.0995 

3 

KCl(cr) 

0.00948 

-0.8002 

4 

MgCl2+12H20(cr) 

0.05052 

1.2364 

5 

Na2SO4*10H2O(cr) 

0.02925 

-12.2171 

AQUEOUS  SOLUTION 

Ionic  strength  11.0759 

Osmotic  coefficient  2.0008 


N  Species 

Moles 

Molality 

Activity 

Act.  coef. 

-G/RT 

1  Na%q) 

0.22059D-02 

0.53403D+00 

0.23599D+00 

0.4419 

0.0000 

2  K+(aq) 

0.11533D-02 

0.27920D+00 

0.15016D-01 

0.0538 

0.0000 

3  Cat2 

0.94653D-02 

0.22914D+01 

0.15724D+01 

0.6862 

0.0000 

4  Mg+2 

0.46407D-02 

0.11235D+01 

0.65163D+00 

0.5800 

0.0000 

5  Cl-(a3) 

0.31425D-01 

0.76077D+01 

0.29938D+02 

3.9352 

0.0000 

6  SO4-2  (aq) 

0.72974D-04 

0.17666D-01 

0.64149D-02 

0.3631 

0.0000 

7  CaS04(aq) 

0.64693D-04 

0.15662D-01 

0.15662D-01 

1.0000 

-0.4410 

8  MgS04(aq) 

0.87333D-07 

0.21142D-04 

0.21142D-04 

1.0000 

5.2851 

9  H?0(1) 

0.22929D+00 

0.65191D+00 

0.0000 

BALANCE 

Total 

Solids 

Solution 

Total  computed 

h2o 

0.555084E+02 

0.552791E+02 

0.229289E+00 

0.555084E+02 

Na+ 

0.486950E+00 

0.484744E+00 

0.220592E-02 

0.486950E+00 

K+ 

0.106300E-01 

0.947671E-02 

0.115329E-02 

0.106300E-01 

Ca+2 

0.953000E-02 

0.000000E+00 

0.953000E-02 

0.953000E-02 

Mg+2 

0.551600E-01 

0.505192E-01 

0.464076E-02 

0.551600E-01 

Cl- 

0.568180E+00 

0.536755E+00 

0.314254E-01 

0.568180E+00 

so4-2 

0.293900E-01 

0.292522E-01 

0.137754E-03 

0.293900E-01 

Number  of  iterations  93 

Temperature 

-55.00°C  (218.15  K) 

SOLID  PHASES 

N  Phase 

Moles 

-G/RT 

1  H20(cr,I) 

53.63903  -0.5004 

2  NaCl*2H20(cr)  0.42817  0.8123 

3  KCl(cr) 

0.01063  -1.4388 

4  CaCl2*6H20 

fcr)  0.00953  5.6614 

5  MgCl2H2H20(cr)  0.05516  -1.5678 

6  Na2SO4n0H2O(cr))  0.02939  -15.0450 

BALANCE 

Total 

Solids 

Solution 

Total  computed 

h2o 

0.555084E+02 

0.555084E+02 

0.0000O0E+00 

0.555084E+02 

Na+ 

0.486950E+00 

0.486950E+00 

0.000000E+00 

0.486950E+00 

K+ 

0.106300E-01 

0.106300E-01 

0.000000E+00 

0.106300E-01 

Ca+2 

0.953000E-02 

0.953000E-02 

0.000000E+00 

0.953000E-02 

Mg+2 

0.551600E-01 

0.551600E-01 

0.000000E+00 

0.551 600E-01 

Cl- 

0.568180E+00 

0.5681 80E+00 

0.000000E+00 

0.5681 80E+00 

so4-2 

0.293900E-01 

0.293900E-01 

0.000000E+00 

0.293900E-01 

Number  of  iterations  88 
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Table  3.  FREZCHEM2  model  output  for  evaporation  of  seawater  at  0°C. 


Temperature  0.00°C  (273.15  K) 

Water  amount  50.00  g 

SOLID  PHASES 

N _ Phase _ Moles  -G/RT 

1  NaCl(cr)  0.28103  3.4448 

2  Na2SO4*10H2O(cr)  0.02021  -5.7260 


AQUEOUS  SOLUTION 

Ionic  strength  8.0709 

Osmotic  coefficient  1.5573 


N 

Species 

Moles 

Molality 

Activity 

Act.  coef. 

-G/RT 

1 

Na+(aq) 

0.16551D+00 

0.35701D+01 

0.28612D+01 

0.8014 

0.0000 

2 

K+(aq, 

0.10630D-01 

0.22929D+00 

0.74428D-01 

0.3246 

0.0000 

3 

Ca+2 

0.82305D-02 

0.17754D+00 

0.16740D+00 

0.9429 

0.0000 

4 

Mg+2 

0.55155D-01 

0.11897D+01 

0.27250D+01 

2.2905 

0.0000 

5 

Cl  (aq) 

0.28715D+00 

0.61941D+01 

0.10952D+02 

1.7681 

0.0000 

6 

S04-2(aq) 

0.78786D-02 

0.16995D+00 

0.10190D-01 

0.0600 

0.0000 

7 

CaS04(aq) 

0.12993D-  02 

0.28027D-01 

0.28027D-01 

1.0000 

-2.7979 

8 

MgS04(aq) 

0.53787D-05 

0.11602D-03 

0.11602D-03 

1.0000 

5.4788 

9 

h2o(1) 

0.25733D+01 

0.72306D+00 

0.0000 

BALANCE 


Total 

Solids 

Solution 

Total  computed 

h2o 

0.277542E+01 

0.202067E+00 

0.257334E+01 

0.277541E+01 

Na+ 

0.486950E+00 

0.321441E+00 

0.165509E+00 

0.486950E+00 

K+ 

0.106300E-01 

O.OOOOOOE+OO 

0.106300E-01 

0.106300E-01 

Ca+2 

0.953000E-02 

0.000000E+00 

0.952987E-02 

0.952987E-02 

Mg+2 

0.551 600E-01 

O.OOOOOOE+OO 

0.551600E-01 

0.551600E-01 

Cl- 

0.568180E+00 

0.281028E+00 

0.287154E+00 

0.568182E+00 

so4-2 

0.293900E-01 

0.202067E-01 

0.91 8331 E-02 

0.293900E-01 

Number  of  iterations  45 


Table  4.  Temperatures  (°C)  of  first  appearance 
of  solid  phases  on  chilling  of  seawater. 

FREZCHEM2 


Solid 

Experiment 

Model a 

model 

Ice 

-1.921b 

-1.924 

-1.921 

Mirabilite 

-8.2C 

-5.90 

-5.87 

Hydrohalite 

-22.9  c 

-22.84 

-22.87 

Sylvite 

-36  c 

-34.25 

-34.30 

MgCl2  *  12H20 

-36 c 

-36.82 

-36.82 

Antarcticite 

-54  c 

-53.64d 

-53.73 

a.  Spencer  et  al.  (1990). 

b.  Fujino  et  al.  (1974). 

c.  Nelson  and  Thompson  (1954). 

d.  Calculated  in  sulfate-free  system. 
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APPENDIX  A:  FORTRAN  LISTING  OF  THE  FREZCHEM2  PROGRAM 


PROGRAM  READWRITE 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CHARACTER  *40  NMC , NEL , NM, NMCOMP , TITLE 

DIMENSION  B AL (10) , NEL (10) , GO ( 40 , 6 ) , G ( 40 ) , GIN (40) ,BALS(10) 
DIMENSION  BALL (10), ACT (20) 

DIMENSION  BC (10,20) ,A(10,40) , EX (40) , NMC (20) /NM(40) , AIN (10, 40) 
DIMENSION  IZC (20) , IZ(20) /NUM(40) /NUMC(20) /NUMIN(40) ,NMM{40) 

COMMON  / COMPNT/  NUR ( 10 ) , IP, NUL ( 10 ) , IPNUL , KCOMP 
COMMON  /NUMBER/  NCAT , NANI ,  NNEI , NMM 
COMMON  /MASSl/  A ,  G , NUM ,  X  ( 4  0  ) 

COMMON  /MASS2 /  NM 
COMMON  /INTGER/  NL,N,N1 

COMMON  /MATRIX/  ACT, UM,  PHI , AH20, IZ, IPP, NMV 
COMMON  /AIN/  AIN , GIN , T , NUMIN , IV 
COMMON  /OUTPUT/  JOPEN ( 10 ) , IOPEN 

10  FORMAT (A40) 

11  FORMAT ( A4  0 / ) 

12  FORMAT (13/) 

13  FORMAT (A15, 13, IX, I2,A40,F12 . 6) 

14  FORMAT ('  THE  CHARGE  OF  1  KG  OF  THE  SOLUTION  IS  EQUAL  TO 
$ , F8 . 5 ) 

15  FORMAT(/13X, 'Temperature  ' , F7 . 2 , '  C  ( ' , F7 . 2 , '  K) ' ) 

1 6  FORMAT  ( '  WATER  AMOUNT  ' , F7 . 2 , '  GRAM ' / 

$  '  - ') 

17  FORMAT ( 2 5X, 'SOLID  PHASES' /IX, 'N  Phase  ', 

$  16X, '  Moles  -G/RT') 

1 8  FORMAT ( 12 , IX , A2  3 , F9 . 5 , 2X , F12 . 4 ) 

19  FORMAT (2 IX, 'Ionic  strength  ',F8.4/ 

$  21X, 'Osmotic  coefficient  ',F7.4) 

20  FORMAT ( IX, 'N' ,2X, 'Species  ' ,4X, 'Moles' ,7X, 

$  'Molality  ', 3X, ' Activity  ','Act.coef.  -G/RT') 

2 1  FORMAT ( 12 , IX, A16 , 3D12 . 5 , 2X, F8 . 4 , F10 . 4 ) 

22  FORMAT (12, 1X,A16,D12 .5, 12X,D12 . 5 , 8X, F10 . 4 , ) 

2  3  FORMAT  ( 2  6X ,  ' BALANCE ' / 

$  13X, ' Total 9X, ' Solids ', 8X, ' Solution ', 5X, ' Total  computed') 

24  FORMAT ( IX , A8 , 4 (2X,E12 . 6) ) 

25  FORMAT ('  Number  of  iterations  ',13) 

26  FORMAT ('  THE  PROGRAM  CANNOT  CALCULATE  THIS  CHEMICAL 
$EQUILIBRIUM ' ) 

OPEN ( 1 , FILE= ' DBASE ' ) 

OPEN ( 9 , FILE= ' INPUT ' ) 

OPEN ( 10 , FILE= ' RESULT ' ) 

XW=55. 50837 

READ (9, 10)  TITLE 
WRITE (10, 11)  TITLE 
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C . .  READING  INDEPENDENT  COMPONENTS  AND  GIVING  THE  BALANCE 
READ (1,12)  KCOMP 

C..  THE  WATER  BALANCE  (1  KG  H20  =  55.50837  MOLES  ) 

BAL ( 1 ) =XW 
NEL(1)='H20' 

NUMIN(l) =1 
NMM(l) =30 
NM ( 1 ) =NEL ( 1 ) 

NUR ( 1 ) =1 
J=1 
J1  =  0 

CHARGE=0 . 

DO  2000  1  =  2 , KCOMP 

READ{1,  13)  NMCOMP,  IIZ,NU 
READ (  9  ,  *  )  AA 
IF  (AA.GT.l.D-20)  THEN 
J=J+1 
BAL  ( J )  =  AA 

C  HARGE = C  HARGE +IIZ*AA 
NEL(J) =NMCOMP 
NM ( J) =NEL { J) 

NMM(J) =NU 
NUMIN { J) =J 
NUR ( J ) =1 

ELSE 

J1=J1+1 
NUL ( J1 ) =1 
END  IF 
2000  CONTINUE 

IF (DABS (CHARGE) .GT. 0.00001)  THEN 

WRITE (*,*)  'THE  SALT  BALANCE  IS  NOT  CORRECT ' 

WRITE (  * , 14 )  CHARGE 

IF ( CHARGE . GT . 0 . ) WRITE (*,*)' ADD  ANIONS  OR  SUBT .  CATIONS  IN 
$ INPUT  7 

IF ( CHARGE . LT . 0 . ) WRITE (  *  ,  *  ) ' ADD  CATIONS  OR  SUBT.  ANIONS  IN 
$ INPUT' 

PAUSE 
GO  TO  815 
END  IF 


IP=J 
IPP=IP 
IPNUL= Jl 

DO  2100  1=1 , IP 

DO  2101  J=l, IP 
AIN(I, J)=0. 
2101  CONTINUE 

AIN (I, I)=l. 

GIN ( I) =0 . 

2100  CONTINUE 
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2201 


2200 


2301 


2300 


2401 


2400 


2501 


CALL  CHOICE  (NCAT , BC , NMC , NUMC , IZC) 
NC=IP+NCAT 

DO  2200  1=1,  NCAT 
I1=I+IP 


DO  2201  J=l, IP 

AIN ( J, II) =BC (J,  I) 

CONTINUE 

NM (II) =NMC ( I ) 

NMM (II) =NUMC ( I ) 

NUMIN(Il) =11 
IZ(I1)=IZC(I) 

CONTINUE 

CALL  CHOICE  (NANI , BC , NMC , NUMC , IZC ) 
NA=NC+NANI 

DO  2300  1=1,  NANI 
I1=I+NC 


DO  2301  J=l, IP 

AIN ( J , II) =BC ( J , I ) 

CONTINUE 

NM  (II)  =NMC  ( I ) 

NMM  (II)  =NUMC  (I) 

NUMIN(Il) =11 
IZ (II) =IZC (I) 

CONTINUE 

CALL  CHOICE  (NNEI , BC , NMC , NUMC , IZC) 
NL=NA+NNEI 

DO  2400  1=1,  NNEI 
I1=I+NA 


DO  2401  J=l, IP 

AIN ( J, II ) =BC ( J, I ) 

CONTINUE 

NM(I1)  =NMC  (I) 

NMM { I 1 ) =NUMC ( I ) 

IF (NMM (II) . EQ . 30 ) IV=I1 
NUMIN(Il) =11 
IZ (II) =IZC (I) 

CONTINUE 

CALL  CHOICE  (NSOL , BC , NMC , NUMC , IZC ) 

DO  2500  1=1,  NSOL 
I1=I+NL 


DO  2501  J=l, IP 

AIN ( J , II) =BC (J, I) 
CONTINUE 
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2500 


2600 

2700 
C.  .  . 

2810 


2820 

2800 


1 

2900 


NM  (II)  =NMC  ( I ) 
NMM(Il) =NUMC (I) 
NUMIN(Il) =11 
CONTINUE 


IP1=IP+1 

L=NCAT+NANI+NNEI 

N=NL+NSOL 

N1=N+1 


DO  2600  J=1,IP 

AIN(J,N1) =BAL (J) 
CONTINUE 


IJJ=0 


DO  2700  1=1,  Nl 
NUM ( I ) =NUMIN ( I ) 
CONTINUE 


READ  TABLE  3 


OPEN (5,  FILE='TABLE3 ' ) 

DO  2800  I=IP+1 ,  N 

READ ( 5  ,  *  ,  END=2  82  0 )  NTl 

IF (NMM { I ) . EQ .NTl)  THEN 
BACKSPACE (5) 

READ ( 5 , * )  NTl,  (GO ( I , K ) , K=1 , 6 ) 

REWIND { 5 ) 

ELSE 

GO  TO  2810 
END  IF 

REWIND { 5 ) 

CONTINUE 

CLOSE  (5) 

CALL  PITZER  (T,0,EX) 

READ ( 9 , * )  TINIT 
READ (9,*)  IPATH 

IF  ( IPATH. EQ.l)  THEN 
READ (9,*)  TFIN 
READ ( 9 , * )  DT 
ELSE 

READ { 9 , * )  WFIN 
READ ( 9 , * )  DW 
END  IF 

T=TINIT 
BALWAT= 1000. 

CONTINUE 

DO  2900  I=IP+1 ,  N 

GIN ( I ) =  PF{T,G0 (I, 1) ,G0 (I, 2) ,G0 {I, 3) , GO (I, 4) , GO (I, 5) , 
$  GO (1,6) ) 

CONTINUE 
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3000 

C.  .  . 


3100 

4 

3200 

3301 

3300 


IF ( I PATH . EQ . 1 )  WRITE (*,*)  ' T= ' ,  T 

PHI  =  1 . 

CALL  PITZER  (Tfl,EX) 

CONTINUE 

IF ( I PATH . EQ . 2 )  WRITE  (*,*)  ' AMOUNT  OF  WATER= 1 ,  BALWAT 

CALL  SOL  (ICK,EX, ISOLU, I PATH) 

TC=T-273 . 15 
WRITE (10 ,15)  TC,T 

IF { I PATH . EQ . 2 )  WRITE (10, 16)  BALWAT 

IF ( ICK. GT .400)  GO  TO  5 

IF ( I OPEN . GT . 0 )  THEN 
WRITE (10, 17) 


DO  3000  1=1 , IOPEN 
J= JOPEN ( I ) 

WRITE (10, 18)  I , NM ( NUM ( J ) ) ,X(J) ,G(J) 
CONTINUE 

END  IF 

IF ( I SOLU . EQ . 0 )  GO  TO  4 


.  START  OF  OUTPUT  . 

WRITE (10,*)  '  AQUEOUS  SOLUTION1 

WRITE ( 10 , 19 )  UM,  PHI 
WRITE (10, 20) 

DO  3100  K=IP1 , NL-1 
EX (K) =EX (K) *XW 
G (K) =G(K) -DLOG(XW) 

AC=EXP (ACT (K) ) 

WRITE ( 10 , 21)  K-IP, NM (K) ,X(K) ,EX(K) , EX (K) *AC , AC , G (K) 
CONTINUE 

WRITE(10, 22)  NL-IP, NM (NL) ,X(NL) ,DEXP(AH20) , G (NL) 

WRITE (10, 23) 

DO  3200  K=1 , IP 
BALL ( K ) =  0 . 

BALS (K) =0 . 

CONTINUE 

DO  3300  1=1, IOPEN 
J= JOPEN { I ) 

DO  3301  K=l, IP 

BALS (K) =BALS (K) +X(J) *AIN (K, NUM ( J) ) 

CONTINUE 

CONTINUE 
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DO  3400  1=1,1? 


DO  3401  J=IP1,NL 

BALL ( I ) =BALL ( I ) +X ( J ) * AIN ( I , J ) 

3401  CONTINUE 

WRITE (10, 24) NEL ( I ) , AIN ( I , Nl ) , BALS(I) , BALL (I) , BALS ( I ) +BALL ( I ) 
3400  CONTINUE 


5  WRITE (10,25)  ICK 

C .  FINISH 


IF ( ICK . GT .400)  WRITE  (10,26) 

IF  ( IPATH . EQ . 1 )  THEN 
T=T-DT 

IF (T . GE. TFIN- .001)  GO  TO  1 
ELSE 

B AL WAT = B ALWAT - DW 

IF ( B ALWAT . GE. WFIN)  THEN 

AIN ( 1 , Nl ) =BALWAT/ 18 .0153 
GO  TO  2 
END  IF 
END  IF 

CLOSE (3) 

815  STOP 
END 


c - C 

c - c 

SUBROUTINE  CHOICE  (NCAT , BK, NMCUR, NUMCUR, IZC) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C..  SUBROUTINE  FOR  CHOOSING  COMPONENTS  FROM  DATABASE  ACCORDING  TO 
C . .  THE  GIVEN  INDEPENDENT  COMPONENTS 

CHARACTER  *40  NM, NMCUR 

DIMENSION  BK ( 10 , 20) , NMCUR (20) , NUMCUR (20) ,IZC(20),B(10) 

COMMON  /COMPNT/  NUR ( 10 ) , IP , NUL ( 1 0 ) , IPNUL , KCOMP 

3  FORMAT (A23 , 12 , IX, 12 , IX, 10F6 . 2 ) 


READ { 1 , * ) NCAT 
J=0 

DO  2  1=1, NCAT 

C  READ (1,3)  NM , NUM , I Z ,  (B( II), 11=1,  KCOMP) 

READ(1,*)  NM , NUM , I Z ,  (B ( II), 11=1,  KCOMP) 

DO  2000  Jl=l, IPNUL 

IF  ( B (NUL ( Jl ) ) . NE . 0 . )  GO  TO  2 
2000  CONTINUE 
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IS  =  0 


DO  3000  Jl=l, IP 

IF  ( DABS (B (NUR ( J1 ) ) )  . GT . ID- 2  0 )  IS=IS+1 
3000  CONTINUE 

IF  (IS.EQ.O)  GO  TO  2 
J=J+1 

DO  4  Jl=l, IP 

BK(J1, J) =B (NUR ( J1 ) ) 

4  CONTINUE 

NMCUR ( J ) =NM 
NUMCUR(J) =NUM 
IZC(J)=IZ 
2  CONTINUE 

NCAT=J 

RETURN 

END 


C - C 

SUBROUTINE  INTACT ( Z1 , Z2 , UM, A, PHIPHI , PHIPRI , PHIIJ, THETA) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C  THIS  SUBROUTINE  CALCULATES  THE  HIGHER-ORDER  ELECTROSTATIC 
C  INTERACTION  TERMS  FOR  PITZER  EQUATIONS. 

DIMENSION  B (0:22) ,  D(0:22),  ZA(2),  XA (2 , 2 ) , XJ (2 , 2 ) , XJPRIM ( 2 , 2 ) 

DIMENSION  AKI (0:20) ,  AKII(0:20) 

DATA  AKI/1. 925154014814667,  -0.060076477753119, 

$  -0.029779077456514,  -0.007299499690937,  0.000388260636404, 

$  0.000636874599598,  0.000036583601823,  -0.000045036975204, 

$  -0.000004537895710,  0.000002937706971,  0.000000396566462, 

$  -0.000000202099617,  -0.000000025267769,  0.000000013522610, 

$  0.000000001229405,  -0.000000000821969,  -0.000000000050847, 

$  0.000000000046333,  0.000000000001943,  -0.000000000002563, 

$  -0.000000000010991/ 

DATA  AKII/0. 628023320520852,  0.462762985338493, 

$  0.150044637187895,  -0.028796057604906,  -0.036552745910311, 

$  -0.001668087945272,  0.006519840398744,  0.001130378079086, 

$  -0.000887171310131,-0.000242107641309,  0.000087294451594, 

$  0.000034682122751,-0.000004583768938,  -0.000003548684306, 

$  -0.000000250453880,  0.000000216991779,  0.000000080779570, 

$  0.000000004558555,-0.000000006944757,  -0.000000002849257, 

$  0.000000000237816/ 

B  ( 21 ) =0 . 

B  ( 22 ) =0 . 

D(21) =0 . 

D ( 22 ) =0  . 

ZA ( 1 ) =Zl 

ZA(2 ) =Z2 

SQ=SQRT(UM) 
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DO  2000  J=1 , 2 


DO  2001  I  =1,2 

XA(J,I)=6.*ZA(I)*ZA(J) *A*SQ 
X=XA( J, I) 

IF  (X.LT.l.)  THEN 

ZZ=4 . *X**0 . 2-2 . 0 
DZ= . 8*X*  * ( - . 8 ) 


DO  2002  K=20 , 0,-1 

B (K) =ZZ*B (K+l) -B (K+2) +AKI (K) 

D  (K) =B (K+l ) +ZZ*D (K+l ) -D (K+2 ) 

2002  CONTINUE 

ELSE 

ZZ=40./9.*X**(-.l)-22./9. 

DZ=-40./90.*X** (-1.1) 

DO  2003  K=20 ,0,-1 

B (K) =ZZ*B (K+l) -B (K+2) +AKII (K) 

D (K) =B (K+l ) +ZZ*D (K+l ) -D (K+2 ) 

2003  CONTINUE 

END  IF 

XJ(J, I)=. 25*X-1. + .5* (B ( 0 ) -B ( 2 ) ) 

XJPRIM (J,I)=.25+.5*DZ* (D ( 0 ) -D ( 2 ) ) 

2001  CONTINUE 

2000  CONTINUE 

ETHETA= (Zl*Z2/4. /UM) * (XJ(1, 2) - . 5*XJ(1, 1) - . 5*XJ(2, 2) ) 
ETHPRI=-ETHETA/UM+ (Zl*Z2/8 . /UM**2) * (XA (1 , 2 ) *XJPRIM ( 1 , 2 ) - 
$  . 5*XA(1, 1) *XJPRIM(1, 1) - . 5*XA (2,2) *XJPRIM(2,2) ) 

PHIPHI=THETA+ETHETA+UM*ETHPRI 
PHI I J=THETA+ETHETA 
PHIPRI=ETHPRI 

RETURN 

END 


C - C 

c - c 

SUBROUTINE  PITZER  (T, IFLAG, EX) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  A0M ( 6 ) ,NMM(40) ,NAQ(20) ,Z(20),EX(20) 

DIMENSION  BPRIME(10, 10) , BPHI (10, 10),C(10,10), PHIIJ(10, 10) , 

$  SUMCAdO)  ,  SUMCAT  (10)  ,SUMAN(10)  ,  SUMZ(10)  ,SUMAC(10)  ,ACT(20)  , 
$  PHIPHI (10, 10) , PHIPRI (10, 10) , XM ( 40 ) , IZ ( 20 ) 

DIMENSION  SUMAA(IO) ,SUMCC(10) , SUMK(10) 

DIMENSION  BET0 (10, 10) , BET1 (10, 10) ,C0 (10, 10) ,  BET2(10,10), 

$  TET (10 , 10 ) ,PSI(10,10,10) ,B(10,10) 

DIMENSION  BETOMdO,  10, 6)  ,BET1M(10, 10,6)  ,  COM  (10, 10, 6)  , 

$  BET2M(10, 10, 6) ,TETM(10, 10,6), PSIM(10, 10,10,6) 

COMMON  /NUMBER/NCAT , NANI , NNEI , NMM 
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COMMON  /MATRIX/  ACT , UM, PHI , AH20 , IZ , IP , NMV 

IF  (IFLAG.EQ.l)  GO  TO  1 
IF  (IFLAG.EQ.2)  GO  TO  2 

ICOL=6 

NN=NCAT+NANI 

NC1=NCAT+1 

DO  2000  1  =  1 , NN 

Z(I)=IZ(I+IP) 

NAQ ( I ) =NMM ( I+IP) 

2000  CONTINUE 

C .  READ  TABLES  1  AND  2 


OPEN (3 ,  FILE= ' TABLE 1 ' ) 


C 


READ  TABLE  1 


DO  3000  1=1, NN 

DO  3001  J=1 , NN 

DO  3002  K=1 , 6 

COM ( I , J, K) =0 . 
BET2M ( I , J, K) =0 . 
3002  CONTINUE 

3001  CONTINUE 

3000  CONTINUE 

READ (  3  ,  *  )  ( A0M ( I ) ,  1=1,6) 


DO  4000  1=1 , NCAT 

DO  4001  J=NC1 , NN 

101  READ ( 3 , FMT= ' (213) '  )  NTl , NT2 

IF {NAQ (I) . EQ . NTl . AND . NAQ ( J ) . EQ . NT2 )  THEN 
READ (3, *)  (BET0M ( I , J, K) ,K=1, 6) 

READ (3, *)  (BETlM ( I , J, K) ,K=1, 6) 

IF ( (NTl . EQ . 3  .OR.  NTl.EQ.4)  .AND.  NT2 . EQ . 12 ) THEN 
READ (3, *)  (BET2M ( I , J, K) ,K=1, 6) 

READ (3, *)  (COM (I, J,K) ,K=1, 6) 

ELSE 

READ ( 3 , * )  (COM { I , J, K) ,K=1, 6) 

END  IF 

ELSE 

IF ( (NTl . EQ . 3  .OR.  NTl.EQ.4)  .AND.  NT2 . EQ . 12 ) THEN 
READ ( 3 , FMT= '  (  /  /  /  )  ' ,END=910) 

ELSE 

READ ( 3 , FMT=  7  (//)  ' ,END=910) 

END  IF 
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GO  TO  101 


4001 

4000 

C.  .  . 


102 


5001 

5000 


105 


6002 

6001 

6000 

103 


END  IF 
CONTINUE 

CONTINUE 

CLOSE (3) 


READ  TABLE  2 


OPEN ( 4 ,  FILE= ' TABLE 2 1 ) 

DO  5000  1  =  1 , NCAT-1 

DO  5001  J=I+1/NCAT 

READ ( 4 , * )  NTl,NT2,NT3 

IF (NAQ (I) .EQ.NTl.AND.NAQ(J) .EQ.NT2 . AND . NT3 . EQ . 0 ) THEN 
BACKSPACE (4) 

READ { 4 ,  *)  NTl,NT2,NT3, (TETM ( I ,  J,  K) ,K=1, 6) 

REWIND ( 4 ) 

ELSE 

GO  TO  102 
END  IF 

CONTINUE 

CONTINUE 

DO  6000  1=1, NCAT-1 

DO  6001  J=I+1 , NCAT 

DO  6002  IJ  =  NC 1 , NN 

READ ( 4 , * )  NTl , NT 2 , NT3 

IF (NAQ (I)  .EQ.  NTl 

$  .AND.  NAQ { J )  .EQ.  NT 2 

$  .AND.  NAQ ( IJ)  .EQ.  NT3 )  THEN 

BACKSPACE (4) 

READ ( 4 , * ) NTl , NT2 , NT3 , ( PSIM ( I , J , I J , K ) , K= 1 , 6 ) 
REWIND (4) 

ELSE 

GO  TO  105 
END  IF 

CONTINUE 

CONTINUE 

CONTINUE 

DO  7000  I=NC1,NN-1 

DO  7001  J=I+1,NN 

READ { 4  ,  *  )  NTl , NT2 , NT3 
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$ 

$ 


IF (NAQ ( I )  .EQ.  NTl 

.AND.  NAQ(J)  .EQ.  NT 2 
.AND.  NT3.EQ.0)  THEN 
BACKSPACE ( 4 ) 

READ (4, *)  NTl , NT2 , NT3 , ( TETM ( I , J , K ) ,K=1, 6) 

REWIND (4) 

ELSE 

GO  TO  103 
END  IF 

7001  CONTINUE 

7000  CONTINUE 

DO  8000  I=NC1 , NN-1 

DO  8001  J=I+1 , NN 

DO  8002  I  A=  1 , NCAT 

106  READ ( 4 , * )  NTl , NT2 , NT3 

IF (NAQ (I)  .EQ.  NTl 

.AND.  NAQ ( J)  .EQ.  NT 2 

.AND.  NAQ(IA)  . EQ .  NT3 )  THEN 
BACKSPACE (4) 

READ ( 4 , * ) NTl , NT2 , NT3 , ( PSIM ( I , J, IA,K) ,K=1, 6) 
REWIND ( 4 ) 

ELSE 

GO  TO  106 
END  IF 


8002  CONTINUE 

8001  CONTINUE 

8000  CONTINUE 
104  CLOSE (4) 

C .  END  OF  READING 

RETURN 

1  CONTINUE 


C . .  CALCULATION  OF  VALUES  OF  THE  PARAMETERES  FOR  CURRENT  TEMPERATURE 
A0  =  PF (T , A0M ( 1 ) , A0M ( 2 ) ,A0M(3) ,A0M(4) ,A0M(5) ,A0M(6) ) 

DO  9000  1=1, NCAT 

DO  9001  J=NC1 , NN 

BET0 (I, J) =  PF (T , BET0M ( I , J, 1) , BET0M ( I , J, 2 ) , 

$BET0M ( I , J, 3 ) , BET0M (I , J, 4) , BET0M ( I , J, 5 ) , BET0M ( I , J, 6 )  ) 

BET1 (I, J) =  PF ( T , BET 1M ( I , J , 1) , BET1M (I , J, 2 ) , 

$BET1M(I, J, 3) , BET1M (I , J, 4) , BET1M ( I , J, 5 ) , BET1M ( I , J, 6 )  ) 

BET2 ( I , J) =  PF (T , BET2M ( I , J, 1 ) , BET2M (I , J, 2 ) , 

$BET2M ( I , J, 3 ) , BET2M ( I , J, 4 ) , BET2M ( I , J, 5 ) , BET2M ( I , J, 6 )  ) 

C0(I,J)=  PF (T, COM (I , J, 1) , COM ( I , J, 2 ) ,C0M(I,J,3) , 


$ 

$ 
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9001 

9000 


10002 

10001 

10000 

2 


11000 


12000 


C0M{I,J,4)  ,C0M(IfJ,5)  fC0M(I,J,6)  ) 

CONTINUE 

CONTINUE 

DO  10000  1=1 , NN-1 

DO  10001  J=I+1,NN 

TET ( I , J) =  PF(T/TETM(I/ J, 1) , TETM ( I , J , 2 ) , TETM ( I ,  J ,  3  ) , 
:  TETM (I, J, 4) , TETM { I ,  J ,  5  ) fTETM(I# J, 6) ) 


DO  10002  IA=1,NN 

PSI ( I ,  J, IA)=PF (T, PSIM ( I , J, IA, 1 ) , PSIM ( I , J, IA, 2 ) , 
PSIM(I, J, IAf 3) , PSIM (I, J, IA, 4) , 
PSIMd,  J,  IA,  5)  ,  PSIM  (I,  J,  IA,  6)  ) 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

CONTINUE 

DO  11000  1=1,  10 


SUMCA ( I ) 

= 

0 

SUMCAT ( I ) 

= 

0 

SUMAN ( I ) 

= 

0 

SUMZ ( I ) 

= 

0 

SUMAC ( I ) 

= 

0 

SUMAA(I) 

= 

0 

SUMCC ( I ) 

= 

0 

SUMK(I) 

= 

0 

CONTINUE 


SMX=0 
ZZ=0  . 
UM=0 . 


DO  12000  I=IP+1 , NMV 
J=I-IP 

XM(J)=EX(I)*55. 50837 
ZZ=ZZ+XM ( J) *DABS ( Z ( J) ) 
UM=UM+XM(J) *Z(J) **2 
SMX=SMX+XM(J) 

CONTINUE 


UM=UM/2 . 

SQ=SQRT{UM) 

ALPHA=2 . *SQ 
ALPHA1=1 . 4*SQ 
ALPHA2=12 . *SQ 

Gl  =  2 . * ( 1 . - { 1 . +ALPHA1 ) *EXP ( -ALPHAl )  ) /ALPHAl *  *  2 
G2=2 . * ( 1 . - ( 1 . +ALPHA2 ) *EXP (-ALPHA2 ) ) /ALPHA2**2 

GPRIl=~2 . * (1 (1. +ALPHAl+ALPHAl**2/2 . ) *EXP ( -ALPHAl ) ) /ALPHA1**2 
GPRI2  =  -2 . * (1. - (1. +ALPHA2+ALPHA2**2/2 . ) *EXP ( -ALPHA2 ) } /ALPHA2**2 
G=2  .  *  ( 1-  ( 1+ALPHA)  *EXP( -ALPHA)  )  /ALP^**2 

GPRIME=-2  .  *  (1 .  -  (1 .  +ALPHA+ ALPHA*  *2/2  .  )  *  EXP  (-ALPHA)  )  /  ALPHA*  *2 


22 


DO  13000  J=1/NCAT 


13001 

13000 


14001 

14000 


15001 
15000 
C.  . 


16002 

16001 


DO  13001  I =NC 1 , NN 

IF  (Z (J) *ABS (Z (I) ) .EQ. 4)  THEN 

BPHI ( J,  I) =BET0 { J, I) +BET1 { J, I) *EXP ( -ALPHAl )  + 

;  BET2 (J, I) *EXP ( -ALPHA2 } 

B (J, I) =BET0 (J, I) +BET1 (J, I) *G1+BET2 (J, I) *G2 
BPRIME(J/ I) =BET1 ( J, I) *GPRIl/UM+BET2 (J, I) *GPRI2/UM 
ELSE 

BPHI (J, I) =BET0 ( J, I) +BET1 ( J, I) *EXP(-ALPHA) 

B ( J7 I) =BET0 ( J, I) +BET1 ( J, I) *G 
BPRIME (J, I) =BET1 (J, I) *GPRIME/UM 
END  IF 

C ( J, I ) =C0 ( J, I) /2 . /SQRT (Z { J) *DABS (Z ( I ) ) ) 

CONTINUE 

CONTINUE 

DO  14000  J=1 , NCAT-1 

DO  14001  1= J+l , NCAT 

CALL  INTACT (Z (J) ,  Z  (I) , UM, A0 , PHIPHI ( J, I) # 

I  PHIPRI ( J, I) , PHIIJ ( J, I) , TET { J , I ) ) 

CONTINUE 

CONTINUE 

DO  15000  J=NC1 , NN-1 

DO  15001  I  =  J  + 1 , NN 

CALL  INTACT {Z (J) # Z (I) ,UM, A0, PHIPHI (J, I) , 

I  PHIPRI (J,  I) , PHIIJ ( J , I) , TET { J , I) ) 

CONTINUE 

CONTINUE 

CALCULATION  OF  SUMMATION  TERMS  FOR  F  AND  PHI. 

SCATON=0 . 

SUBSUM=0 . 

SANON=0 . 

SUMCAF  =  0 . 

SUMANF  =  0 . 

DO  16000  J=l, NCAT-1 

DO  16001  Jl= J+l , NCAT 


DO  16002  I=NC1 , NN 

SUBSUM=SUBSUM+PSI ( J, Jl , I ) *XM ( I ) 

CONTINUE 

SCATON=SCATON+ (SUBSUM+PHIPHI ( J, Jl) ) *XM ( J ) *XM( Jl) 
SUMCAF = SUMCAF + PH I PR I {J, Jl) *XM ( J) *XM(J1) 

SUBSUM=0 . 

CONTINUE 
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16000  CONTINUE 


SUBSUM-0 . 

DO  17000  J=NC1,NN-1 

DO  17001  Jl= J+l , NN 

DO  17002  1=1 , NCAT 

SUBSUM=SUBSUM+PSI (J, Jl , I) *XM(I) 

17002  CONTINUE 

SANON=SANON+ (SUBSUM+PHIPHI ( J, Jl) ) *XM( J) *XM ( Jl) 
SUMANF = SUMANF+ PHI PRI (J, Jl) *XM{ J) *XM( Jl) 
SUBSUM=0 . 

17001  CONTINUE 

17000  CONTINUE 


SUMB=0 . 

SUMPHI=0 . 

DO  18000  J=1 , NCAT 

DO  18001  I=NC1 , NN 

SUMB=SUMB+XM ( J) *XM(I) *BPRIME(J, I) 

SUMPHI=SUMPHI+XM ( J) *XM(I) * (BPHI ( J, I) +ZZ*C ( J, I) ) 
18001  CONTINUE 

18000  CONTINUE 

F=-A0* (SQ/ (1 . +1 . 2*SQ) +2 . *DLOG ( 1 . +1.2*SQ)/1.2)+ 


$  SUMB+SUMCAF+ SUMANF 

PHI  =  1. +2 . / SMX* ( -A0*UM*  *1.5/ (1 . +1 .2*SQ) +SUMPHI  +  SCATON+SANON) 
AH20=-PHI*SMX/55 .50837 

C .  CALCULATION  OF  TERMS  FOR  ACTIVITY  COEFFICIENTS (GAMMA)  ...  . 


SUM=0 . 


DO  19000  J=l,  NCAT-1 

DO  19001  J1=J+1#  NCAT 

DO  19002  I=NC1 , NN 

PSI (Jl, J, I) =PSI (J, Jl, I) 
19002  CONTINUE 

PHIIJ( Jl, J) =PHIIJ(J, Jl) 

19001  CONTINUE 

19000  CONTINUE 

DO  20000  I=NC1 ,  NN-1 

DO  20001  11=1+1,  NN 

DO  20002  J=l,  NCAT 

PSI (II, I, J) =PSI (I, II, J) 
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20002 


CONTINUE 


20001 

20000 

21001 

21000 


22001 

10 

22000 


23002 

23001 

23000 


24001 

24000 

25000 


PHIIJ(I1, I) =PHIIJ(I, II) 

CONTINUE 

CONTINUE 

DO  21000  J=l,  NCAT 

DO  21001  I=NC1 , NN 

SUMCA(J) =  SUMC A ( J ) +XM(I) * (2.*(B(J,I)) +  ZZ*C ( J, I) ) 
CONTINUE 

CONTINUE 

DO  22000  J=l,  NCAT 

DO  10  Jl=l , NCAT 

IF  (J.EQ.Jl)  GO  TO  10 

DO  22001  I=NCl,  NN 

SUM=SUM+XM ( I ) *  PSI < J, Jl, I) 

CONTINUE 

SUMCAT ( J) = SUMC AT (J) +XM( Jl) * (SUM+2 . *PHIIJ(J, Jl) ) 
SUM=0 . 

CONTINUE 

CONTINUE 

DO  23000  J=l,  NCAT 

DO  23001  Jl=NCl , NN- 1 

DO  23002  I=Jl+l,  NN 

SUMAN ( J ) =SUMAN ( J) +XM( Jl) *XM(I) *PSI (Jl, I, J) 
CONTINUE 

CONTINUE 

CONTINUE 


SUM=0 . 


DO  24000  J=1 , NCAT 

DO  24001  I=NC1 , NN 

SUM=  SUM+XM { J ) *  XM ( I ) *C(J, I) 

CONTINUE 

CONTINUE 

DO  25000  J= 1 , NCAT 

SUMZ ( J) =SUM*DABS ( Z { J) ) 

CONTINUE 

DO  26000  J=1 , NCAT 

ACT ( J+IP) =Z ( J) *  * 2  * F + SUMCA ( J ) + SUMCAT ( J) +SUMAN ( J) +SUMZ ( J) 
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26000  CONTINUE 


27001 

27000 


28001 

20 

28000 


29002 

29001 

29000 


30001 

30000 

31000 

32000 


SUM=0 . 


DO  27000  I=NC1,  NN 

DO  27001  J=l,  NCAT 

SUMAC (I) =  SUMAC (I) +XM( J) *(2.*(B(J,I)) +  ZZ*C (J, I) ) 
CONTINUE 

CONTINUE 

DO  28000  I=NC1 ,  NN 

DO  20  Il=NCl , NN 

IF  (I ,EQ. II)  GO  TO  20 

DO  28001  J=l,  NCAT 

SUM=SUM+XM(J) *PSI (I, 11/ J) 

CONTINUE 

SUMAA ( I ) =SUMAA { I ) +XM ( II ) * (SUM+2 . *PHIIJ ( I, II) ) 
SUM=0  . 

CONTINUE 

CONTINUE 

DO  29000  I=NC1 ,  NN 

DO  29001  J=1 , NCAT- 1 

DO  29002  Jl=J+l,  NCAT 

SUMCC(I) =SUMCC (I) +XM(J) *XM(J1) *PSI ( J,  Jl ,  I) 
CONTINUE 

CONTINUE 

CONTINUE 


SUM=  0 . 


DO  30000  J=1 , NCAT 

DO  30001  I=NC1,NN 

SUM=SUM+XM(J) *XM ( I ) *C ( J, I ) 

CONTINUE 

CONTINUE 

DO  31000  I=NC1 , NN 

SUMK { I ) =SUM*ABS { Z ( I ) ) 

CONTINUE 

DO  32000  J=NC1 , NN 

ACT ( J+IP) =Z (J) **2  *F+SUMAC ( J) +SUMAA{J) +SUMCC(J) +SUMK(J) 
CONTINUE 

DO  33000  1=1 , NMV-IP 

IF (ACT ( I+IP ) „ GT . ALOG (10000 . ) )  ACT ( I+IP) =ALOG ( 10000 . ) 
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33000  CONTINUE 


RETURN 

910  WRITE ( * , 1910)  NAQ(I) ,NAQ(J) 

1910  FORMAT (2 (IX, 13) , /, ' ! !  TABLE  1  IS  NOT  COMPLETE  !!'/ 

$ '  ! !  OR  WRONG  COMPONENT  WAS  ENTERED  ! ! ' ) 

STOP 

END 

C - C 

c - c 

SUBROUTINE  SIMPL  { IPl , M, IOPl , IP , Nl , K, IR) 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z) 

COMMON  /MASSl /  A, G, NUM, X { 4 0 ) 

COMMON  /ENTER/  ISIM, ISIMIN 
DIMENSION  A  ( 1 0 , 4  0  )  ,  G  (  4  0  )  , NUM (40) 

C  IPl  -  BEGINING  J#  M  -  FINISH  J,  IOPl-  BEGINING  I,  IP 

C . LOOKING  FOR  A  NEW  PHASE  FOR  INCLUDING  INTO  BASIS... 

ISIM=0 

1111  DEL=0 . DO 
IK=0 
K  =  0 


DO  1001  J=IP1/M 
P=0  . 


DO  1002  1=1,1? 

P  =  P+G { I ) *A{I, J) 

1002  CONTINUE 

DELTA=G ( J) -P 

IF ( DELTA. GE. DEL)  GO  TO  1001 

DEL=DELTA 

K=J 

ISIM=1 

1001  CONTINUE 

IF(K.EQ.O)  GO  TO  1000 

C . LOOKING  FOR  PLACE  IN  THE  BASIS  TO  BE  SUBSTITUTED 

BMIN=lD+2  0 

DO  1005  I=IOPl , IP 

I F  ( A  ( I ,  K ) . LE . 0 . )  GO  TO  1005 
BTEK=A ( I , Nl ) /A ( I , K) 

IF (BTEK . GE . BMIN)  GO  TO  1005 

BMIN=BTEK 

IR=I 

1005  CONTINUE 

DO  1007  1=1, IP 

A  ( I , IR) =A ( I , K) 

1007  CONTINUE 


-FINISH  I 
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NUM(IR)  =NUM(K) 

G(IR)=G(K) 

CALL  GGl ( IR, IP+1 , IP , Nl ) 

C..  IF ( ISIMIN . EQ . 0 )  GO  TO  1111 

1000  CONTINUE 

RETURN 

END 

C— - C 

C— - c 

SUBROUTINE  GGl ( ICOL , IPl , IPS , Nl ) 

DOUBLE  PRECISION  A,G,X 
DIMENSION  A ( 10 , 40 ) , G (40 ) , NUM (40) 

COMMON  /MASS1/  A, G, NUM, X ( 40 ) 

DO  1  I  =  IP1 , Nl 

A ( ICOL , I ) =A ( ICOL , I ) / A ( ICOL , ICOL ) 

1  CONTINUE 

A  ( ICOL, ICOL) =1. 

DO  2  J=1,IPS 

IF ( J . EQ . ICOL)  GO  TO  2 

DO  3  I=IP1,N1 

A( J, I) =A(J, I) -A (ICOL, I ) *A ( J, ICOL) 
3  CONTINUE 

A ( J, ICOL) =0. 

2  CONTINUE 

RETURN 

END 


C- 


■C 


c 


■c 


SUBROUTINE  GG ( INI , NP , NPl ) 

DOUBLE  PRECISION  Pi, AM, AT 
C OMMON / C OMGG /  Pi (2  0, 21) 

DO  201  I=INI , NP 
AM=P1 (I, I) 

11  =  1 
IJ=I+1 

DO  202  K=I J , NP 

IF (ABS (AM) . GE . DABS (Pi (K, I ) ) )  GO  TO  202 
AM=P1 (K, I) 

I1=K 

202  CONTINUE 
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DO  203  K=I , NPl 

AT=P1 (Il,K) 

PI (II , K) =P1 (I,K) 

Pi ( I ,  K) =AT/AM 
203  CONTINUE 

DO  241  J=INI , NP 
AM=P1 (J,  I) 

IF(J.EQ.I)  GO  TO  241 

DO  204  K=I,NP1 

PI ( J#K) =P1 ( J#K) -PI (I,K) *AM 


204 

CONTINUE 

241 

CONTINUE 

201 

CONTINUE 

RETURN 

END 

C - C 

C - C 

C . .  THE  ROUTINE  FOR  CALCULATION  OF  CHEMICAL  EQUILIBRIA 

C..  IN  WATER-SALT  SYSTEMS  AT  THE  TEMPERATURE  RANGE  25-  -60  C. 

DOUBLE  PRECISION  FUNCTION  PF ( T , Al , A2 , A6 , A9 , A3 , A4 ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PF=Al+A2*T+A6*T**2+A9*T**3+A3/T+A4*LOG(T) 

RETURN 

END 

C - C 


c - c 

SUBROUTINE  SOL  ( ICK, EX, ISOLU , I PATH) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CHARACTER  *40  NM 

DIMENSION  G(40) ,  X  ( 4  0 ) , EX (20) /EXTMP(20) ,A(10,40) , ACT (20) /NM(40) 
DIMENSION  JOPEN(IO) , JCLOS(IO) , SU ( 2 0 ) , Pi ( 2 0 , 2 1 ) , AIN (10, 40) 
DIMENSION  NUMIN ( 40 ) /NUM(40) #GIN(40) , ACTT(20) , IZ(20) 

COMMON  /NUMBER/NCAT  ,  NANI , NNEI , NMM ( 40 ) 

COMMON  /MASSl /  A,G,NUM,X 
COMMON  /MASS2 /  NM 
COMMON  /INTGER/  NL, N, Nl 

COMMON  /MATRIX/  ACT , UM, PHI , AH20, IZ , IP , NMV 

COMMON  /AIN/  AIN, GIN, T, NUMIN, IV 

COMMON  /COMGG/  PI 

COMMON  /ENTER/  ISIM, ISIMIN 

COMMON  /OUTPUT/ JOPEN, IOPEN 

DO  1000  I  =  1,  20 
ACT (I)  =  0 . DO 
1000  CONTINUE 
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ISIMIN=1 

IS0LU=1 

ICK=0 

ICE=0 

XW=55. 50837 

DO  2000  1  =  1,  10 

DO  2001  J  =  1,  40 

A ( I , J ) = AIN ( X , J ) 

2001  CONTINUE 

2000  CONTINUE 

DO  2100  1=1,  40 

NUM(I) =NUMIN(I) 

G ( I ) =GIN ( I ) 

2100  CONTINUE 

EPS= . 001 
IP1=IP+1 
NUM(Nl) =N1 
NM (Nl) = ' DEL ' 

PHI=1. 

IP1=IP+1 

NL1=NL+1 

L=NL-IP 

NMV=NL-1 

N1=N+1 

DO  2200  I=IP1 , NMV 

G ( I ) =G ( I ) +DLOG (XW) 

2200  CONTINUE 

DO  2300  1=1, IP 

DO  2301  J=IP1,NL 

IF (NMM(I) .EQ.NMMf J) )  THEN 
NUM ( I ) =J 
NUMIN(I) =J 
G ( I ) =G ( J) 

END  IF 

2301  CONTINUE 


IF(NMM(I) .EQ.30)  THEN 
IW=I 

X  ( IV)  =A  ( I ,  Nl ) 

END  IF 

2300  CONTINUE 


SUM=0 .DO 


DO  2400  1=1, IP 
NI=NUM ( I) 

IF  (I.NE.IW)  THEN 
X (Nl) =A (I ,N1) 

EX (Nl ) =A ( I , Nl ) /X(IV) 

G  ( I ) =G (NUM ( I ) ) +DLOG(EX(NUM(I) ) ) 
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2400 


2500 


2601 


2600 


75 

76 


2700 


2800 


2901 


END  IF 
CONTINUE 

CALL  PITZER  (T,2,EX) 

G  ( IW)  =AH20 

DO  2500  1=1,1? 

NI=NUM ( I ) 

IF  (I.NE.IW)  THEN 

G (I ) =G (NI ) +DLOG ( EX ( NI ) ) +ACT ( NI ) 

END  IF 
CONTINUE 

DO  2600  J=IP1 , NMV 
SU(J)=-G(J) 

DO  2601  K=1 , IP 

SU(J)=SU(J)+A(K, J) *G(K) 

CONTINUE 

EX ( J ) =DEXP ( SU ( J ) -ACT ( J ) ) 

X(J) =X(IV) *EX(J) 

CONTINUE 

CALL  SIMPL  (NLl , NLl , 1 , IP, Nl , IK, IR) 

DELL=1 . 

CONTINUE 

ICLOS=0 

IOPEN=0 

DO  2700  1=1, IP 

IF  (NUM ( I ) . GE . NLl . AND . NUM ( I ) . LE . N)  THEN 
IOPEN=IOPEN+l 
JOPEN ( IOPEN) =1 
ELSE 

ICLOS=ICLOS+l 
JCLOS { ICLOS ) =1 
END  IF 
CONTINUE 

IF  ( IOPEN. GT.0)  THEN 

IF  ( IOPEN. EQ.l. AND. NMM (NUM ( JOPEN (1) ) ) .EQ.31)  THEN 
ICE=1 

SU  ( IV)  =-G  ( IV) 

DO  2800  K=1 , IP 

SU  ( IV)  =SU  { IV)  +A  (K,  IV)  *G(K) 

CONTINUE 

XMN=1 . 

DO  2900  K=1 , 10 

XMN=XMN* 1 . 5 

DO  2901  J=IP1 , NMV 

EXTMP ( J ) =EX { J) *XMN 
SUMS=SUM*XMN 
CONTINUE 
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2900 

10 


3000 


3201 

3200 


3300 

11 

3400 


C 


3500 


CALL  PITZER  (T,2,EXTMP) 

IF  (AH20.LE.SU (IV) )  GO  TO  10 
CONTINUE 

X(IV)=X(IV) /XMN* 1 . 2 

DO  3000  1=1, ICLOS 
J=JCLOS (I) 

NJ=NUM(J) 

EX(NJ) =EXTMP(NJ) 11. 2 
G ( J) =G (NJ) +DLOG (EX (NJ) ) +ACT (NJ) 

CONTINUE 

ELSE 

DO  3200  1=1, IP 

DO  3201  J=IPl,  NL 

IF (NUMIN ( I ) . EQ . J)  THEN 

X(J)=X(J) -.95*A(IR,N1) *AIN(I, IK) 
IF (X ( J) . LE . 0 . ) X ( J) =1 . E-5 
END  IF 
CONTINUE 

CONTINUE 

IF ( ICE . EQ . 1 )  THEN 
SU ( IV) =-G ( IV) 


DO  3300  K=1 , IP 

SU(IV)  =SU  (IV)  +A(K,  IV)  *G  (K) 
CONTINUE 


111=0 
END  IF 

DO  3400  J=IPl , NMV 

EX  ( J)  =X  ( J)  /X  ( IV) 

CONTINUE 

111=111+1 

CALL  PITZER  (T,2,EX) 

IF ( ICE . EQ . 1 . AND . IPATH . EQ . 1 )  THEN 
IF (AH20.GT.SU (IV) ) THEN 
X(IV)=X(IV) /1.2 
GO  TO  11 
ELSE 

IF ( III . GT . 1 ) X ( IV) =X { IV) *1.1 
END  IF 
END  IF 

DO  3500  1=1, ICLOS 
J= JCLOS ( I ) 

NJ=NUM ( J) 

IF (NJ . LT . NL)  G ( J) =G(NJ) +DLOG(EX(NJ) ) +ACT(NJ) 
CONTINUE 
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3601 


3600 


5 

C.  . 


3700 


3801 


3800 


DO  3600  J=IP1 , NMV 

SU(J) =-G ( J) 

DO  3601  K=1 , IP 

SU ( J) =SU { J) +A (K, J) *G(K) 

CONTINUE 

IF ( J .NE . IV)  EX ( J) =DEXP (SU ( J) -ACT ( J) ) 
X(J) =EX ( J) *X(IV) 

CONTINUE 

END  IF 
END  IF 

IN=ICLOS+L 

IN1=IN+1 

IF ( ICK . GT . 400 )  RETURN 
ICK=ICK+1 

EXAMINE  FOR  AQUEOUS  SOLUTION  PRESENCE . 

IF ( ICLOS . LE . 1 )  THEN 
ISOLU=0 

DO  3700  I=IP1 ,NL 
X(I)=0. 

CONTINUE 

GO  TO  400 
END  IF 


ISL  =  0 


DO  3800  J=IP1,NL 
SU ( J) =G ( J) 

DO  3801  K=1,IP 

SU(J) =SU(J) -A  (K,  J) *G ( K) 
CONTINUE 

EX ( J) =X(J) /X(IV) 

ACTT ( J) =ACT ( J) 

IF ( IS . EQ . 1 ) ACTT ( J) =0 . 
CONTINUE 


IS=IS+1 

PHITT=PHI 

CALL  PITZER  (T,2,EX) 

PHI= (PHITT+PHI) /2 
SUM=0 . 

DO  3900  J=IP1 , NMV 

ACT ( J) = (ACT { J) +ACTT (J) ) /2 . 
JI=J-IP 

P1(JI, JI)=1./X(J) 
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3900 


4001 

4000 


4100 


4200 

C.  . 
400 


3101 


Pl(JI,INl)=DLOG(X(J) ) +ACT ( J ) - DLOG ( X ( IV ) )+SU(J) 
PI ( JI , L) =-l . /X (IV) 

Pi (L, JI) =-l. /X(IV) *  PHI 
SUM=SUM+EX(J) 

CONTINUE 

PI <L,L) =SUH*PHI/X ( IV) 

PI (L, INI) =SU(NL) -SUM* PHI 

DO  4000  1=1 , ICLOS 
K=L+I 

II=JCLOS (I) 

Pi  (K,  K)  =0  . 

Pi (K, INI ) =A( II , N1 ) 

DO  4001  Js=IPl,NL 

PI (K, INI ) =P1 (K, INI) -A (II, J) *X(J) 

IJ=J-IP 

PI (K, IJ) =-A ( II , J) 

P1(IJ,K)=-A(II,J) 

CONTINUE 

CONTINUE 

CALL  GG  (1, IN, INI) 

DO  4100  1=1, ICLOS 
K=L+I 

II=JCLOS (I) 

YLI=G(II) -Pi (K, INI) /DELL 
XBS=DABS ( YLI-G (II)  ) 

IF (XBS.LE.EPS)  ISL=ISL+1 
G ( II) =YLI 
CONTINUE 

DO  4200  1=1, L 
J=I+IP 

YLI=X ( J) -PI ( I , INI ) /DELL 
IF (YLI . LE . 1 . D-10 )  YLI=1 . D-6 
XBS=DABS ( (YLI-X(J) ) /X(J) ) 

IF (XBS.LE.EPS)  ISL=ISL+1 
X ( J) =YLI 
CONTINUE 

IF (X (NL ) . LE . 1 . D-3 )  X (NL) =1 . D-3 
IF(ISL.LT.IN)  GO  TO  5 

DELETING  OF  FIXED  COMPOSITION  PHASES 


XX=0. 

IMIN=0 

DO  3100  1=1 , IOPEN 
JO=JOPEN ( I ) 
SUM=0 . 


DO  3101  J=IP1 ,  NL 

SUM=SUM+A( JO, J) *X(J) 
CONTINUE 
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X(JO) =A(JO,Nl) -SUM 
IF(X(JO) .GT.XXJGO  TO  3100 
XX=X(JO) 

IMIN=JO 

3100  CONTINUE 

IF  (ISOLU.EQ.0)  RETURN 
IF  (IMIN.EQ.O)  GO  TO  1025 
NUM(IMIN) =N1 
GO  TO  76 

1025  CALL  SIMPL (NL1 , N, 1 , IP, Nl , IK, IR) 

IF (ISIM.EQ. 1)  THEN 
DELL=DELL+ . 5 
GO  TO  75 
END  IF 

7000  CONTINUE 

RETURN 

END 
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APPENDIX  B:  DATA  FILES  FOR  PROGRAM  FREZCHEM2 

File  DATABASE 

9  -  number  of  independent  components 

H20  0  30  -  the  name  of  the  independent  component,  its  charge  and  its 

Na+  1  1  code  (  the  codes  in  this  model  are  the  same  as 

K+  1  2  in  the  FREZCHEM  model  ) 

Ca2+  2  3 

Mg2+  2  4 

Cl-  -1  11 

S04-2  -2  12 

C03-2  -2  15 

H+  16 


5  -number  of  cations 

H20 

Na  + 

K+ 

Ca2+ 

Mg 2+  Cl- 

S04 

-2  C03 

-2  H+ 

Na+ {aq) 

1 

1 

.00 

1.00 

.00 

.  00 

.00 

.00 

.00 

.00 

.00 

K+ (aq) 

2 

1 

.00 

.00 

1.00 

.00 

.00 

.00 

.00 

.  00 

.00 

Ca2+ (aq) 

3 

2 

.00 

.00 

.00 

1.00 

.00 

.00 

.00 

.  00 

.00 

Mg2  + (aq) 

4 

2 

.00 

.00 

.  00 

.00 

1.00 

.  00 

.00 

.00 

.00 

H+ (aq) 

5  -number  of  anions 

5 

1 

.00 

.00 

.00 
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